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Abstract 

We apply a recently developed perturbative formalism which describes the evolution under their 
self-gravity of particles displaced from a perfect lattice to quantify precisely, up to shell crossing, 
the effects of discreteness in dissipationless cosmological A^-body simulations. We give simple ex- 
pressions, explicitly dependent on the particle density, for the evolution of power in each mode as 
a function of red-shift. For typical starting red-shifts the effect of finite particle number is to slow 
down slightly the growth of power compared to that in the fluid limit (e.g. by about ten percent at 
half the Nyquist frequency), and to induce also dispersion in the growth as a function of direction at 
a comparable level. In the limit that the initial red-shift tends to infinity, at fixed particle density, 
the evolution in fact diverges from that in the fluid limit (described by the Zeldovich approximation). 
Contrary to widely held belief, this means that a simulation started at a red-shift much higher than 
the red-shift of shell crossing actually gives a worse, rather than a better, result. We also study 
how these effects are modified when there is a small-scale regularization of the gravitational force. 
We show that such a smoothing may reduce the anisotropy of the discreteness effects, but it then 
increases their average effect. This behaviour illustrates the fact that the discreteness effects de- 
scribed here are distinct from those usually considered in this context, due to two-body collisions. 
Indeed the characteristic time for divergence from the coUisionless limit is proportional to N^^^, 
rather than N/ log A*' in the latter case. 

PACS numbers: 98.80.-k, 05.70.-a, 02.50.-r, 05.40.-a 



I. INTRODUCTION 



Cosmological A'^-body simulations have become an es- 
sential instrument in the attempt to understand the ori- 
gin of large scale structure in the universe within the 
framework of current cosmological models (see, e.g., [1- 
3]). In this context the goal is to use such simulations to 
recover the clustering of dark matter, which is described 
by a set of coupled Vlasov-Boltzmann equations corre- 
sponding to an appropriately taken N oo limit. The 
problem of discreteness is that of the relation between 
these finite N simulations and this continuum limit. It 
is a problem which is becoming rapidly of greater im- 
portance because of the need for great precision in the 
predictions furnished by simulations, as the observations 
constraining them will continue to improve in the coming 
years (see e.g. [4] for a discussion of precision require- 
ments for future weak lensing surveys). 

This article is the second in a series in which we ad- 
dress this issue in a quantitative manner. In the first 
paper [5] we have considered only the initial conditions 
(IC) for such simulations, studying in detail the correla- 
tion properties of the point distributions generated with 
the standard algorithm used to produce them. We have 



identified the limit in which the continuum IC are recov- 
ered, and quantified precisely the corrections to this limit 
which appear at finite particle number. We have noted 
that there is a subtlety about how this limit is taken: the 
exact theoretical correlation properties of the IC in real 
and reciprocal space are recovered only when the limit 
A'^ — s- oo is taken before the limit that the initial redshift 
•Zinit ^ oo. If, on the other hand, Zinit is increased arbi- 
trarily at fixed particle density, the real space correlation 
properties (e.g. mass variance in spheres) become domi- 
nated by those of the underlying (unperturbed) particle 
distribution, and thus differ completely from those of the 
theoretical IC. As the analysis of [5] is limited to the 
IC, no conclusion can be drawn about the relevance of 
these behaviours in the evolved simulations. One of our 
results here is that this non-interchangeability of these 
limits also manifests itself in the perturbative approxi- 
mation to the evolution we use. As explained briefly in 
the conclusion of [5], this behavior can be understood 
as a manifestation of the fact that the Vlasov limit of 
an A'^-body system with long-range interactions is valid 
for sufficiently short times [6, 7]. More specifically, this 
means that taking Zjnit — » cxd, at fixed particle density, 
the evolution of the particle system diverges from that of 
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the Vlasov-Poisson limit. Practically this implies that in- 
creasing the starting red-shift Zi^it of a simiilation, keep- 
ing all else fixed, the results get worse, in the sense that 
they deviate further from the desired physical behaviour. 
This finding corrects a widespread belief (see, e.g., [8-11]) 
that the opposite is true, as the only envisaged error has 
been that due to non-linear corrections in the fluid limit 
(which vanish when 2;i„it — > c»). 

In this paper we pursue the study of [5] of discrete- 
ness effects, considering the early times dynamics of these 
simulations. By "early times" we mean the regime of va- 
lidity of a perturbativc treatment of the evolution of the 
full A''-body problem which has been reported in detail 
in [12, 13]. This corresponds to a time when the rel- 
ative displacement of particles at adjacent lattice sites 
(which is always small compared to unity in the IC) be- 
come of order the lattice spacing, i.e., when pairs of par- 
ticles first approach one another. In fluid language it 
corresponds approximately to "shell crossing" . The ba- 
sis of the approach is a simple perturbative expansion, 
in the relative displacements of pairs of particles, of the 
force on each particle. This is in fact a standard method 
used to analyze phonons in a crystal, well known in solid 
state physics (see e.g. [14]). It leads to a 3Nx3N diag- 
onalisation problem, which can be simplified using the 
translational symmetry of the lattice (formulated in the 
Bloch theorem): the eigenmodes of the displacement field 
are simply plane waves, with orientations and eigenval- 
\ies which can be determined mimerically by diagonal- 
ising N 3x3 matrices. Analysing this spectrmn in the 
continuum limit of infinite particle density, one recovers 
the expected fluid behaviour in this regime. The latter 
is that derived through a perturbative treatment of the 
fluid equations in the Lagrangian formalism [15], which 
gives asymptotically the Zeldovich approximation. We 
will refer to the discrete perturbative treatment here as 
"particle linear theory" (PLT) and to its fluid limit as 
"fluid linear theory" (FLT). By comparing evolution de- 
scribed by the former with the latter, we determine the 
discreteness effects (in the domain of validity of PLT). 

That this PLT formalism can be used to this end has 
already been clearly demonstrated in [13]. Indeed in 
this paper we have directly compared A'^-body simula- 
tions with both PLT and FLT and shown that the for- 
mer reproduces better than the latter the real evolution. 
Further we have shown that PLT reproduces the full evo- 
lution extremely well at early times, and breaks down 
globally at roughly the same time as FLT, when the av- 
erage relative displacement of particles becomes of order 
the interparticle distance. 

In the present paper we apply this formalism in de- 
tail to quantifying discreteness effects in cosmological A''- 
body simulations. To our knowledge this is the first work 
in which any such effect in simulations has been quanti- 
fied using analytic methods. Discreteness effects have, on 
the other hand, been investigated, mostly numerically, by 
various authors (see e.g. [10, 16-26]). While our results 
are limited to the effects which are at play up to shell 



crossing, they provide a complete and exhaustive under- 
standing of this regime. As we discuss in our conclusions 
the physical insights gained should also be useful in at- 
tempting to understand the effects of discreteness better 
beyond this regime. Our results also provide an analyt- 
ical description for some specific effects of discreteness 
which have been observed numerically, notably the ef- 
fects of discreteness in breaking isotropy demonstrated 
by a numerical experiment in [19], and the generation of 
structure at small scales noted in hot/ warm dark matter 
simulations in [25, 26]. 



The paper is organized as follows. In the next section 
we summarize the necessary formalism which has been 
developed and studied in detail in [12, 13]. This section 
is divided into three parts: the PLT formalism in full gen- 
erality, the derivation of the fluid limit from it, and the 
specific case of evolution in an Einstein de Sitter (EdS) 
universe from the IC applied in cosmological simulations 
(using the Zeldovich approximation). For the latter case 
(i.e. an EdS universe) all of the cosmology dependent 
part of the calculations can be done analytically. It is in 
fact the appropriate case for almost all cosmological ap- 
plications, since we are treating in practice the evolution 
in an epoch (before shell crossing) which corresponds to 
a rod-shift range in which the universe, in current cosmo- 
logical models, is very close to EdS. The generalization 
to any other cosmology (e.g. with a cosmological con- 
stant) is, however, straightforward. In Sect. Ill we then 
apply these results to quantify more precisely the effects 
in cosmological simulations. Specifically we calculate a 
simple function quantifying the modification of the av- 
erage amplitude of all modes of the displacement field 
at given wavenumber, then a similar one for the modes 
of the density fiuctuations and finally one quantifying 
the anisotropy induced in the evolution by discreteness. 
In the following section we study the case that there is a 
smooth regularization of the Newtonian potential around 
the origin. In particular we consider the modification 
induced in the quantities defined and calculated in the 
previous section. In Sect. V we consider the parametric 
and limiting behaviours of the discreteness effects which 
we have quantified. In particular we consider more ex- 
plicitly how the effects depend on the number of particles 
and the recovery of the continuum limit, as well as the 
limit in which the initial amplitude goes to zero (i.e. the 
initial red-shift goes to infinity). In the final section we 
summarize our findings and discuss some other points. 
In particular we discuss the possible use of our results to 
correct cosmological simulations for these systematic er- 
rors due to discreteness at shell crossing. We also discuss 
briefly the relevance of our results to the more general 
problem of quantifying discreteness errors in the fully 
non-linear regime of these simulations. 
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II. PLT EVOLUTION OF A PERTURBED 
LATTICE 

In the first part of this section we present the essential 
elements of the PLT treatment of the early time evolution 
of a perturbed lattice, summarizing a much more detailed 
discussion which can be foimd in [12, 13]. In the next 
subsection we discuss the derivation of the fluid limit of 
the evolution described by PLT. In the last part we then 
consider the specific case of evolution from a lattice as 
perturbed initially in cosmological simulations. 



A. Summary of PLT formalism 

The equation of motion of the particles in a dissipa- 
tionless cosmological A''-body system is (see, e.g., [1, 2]) 



For gravity we have 



Xi + 2H{t)±^ 



in 



Gmj{yii — Xj) 



(1) 



Here dots denote derivatives with respect to time t, Xj 
is the comoving position of the zth particle, of mass m^, 
related to the physical coordinate by = a(t)xj, where 
a{t) is the scale factor of the background cosmology with 
Hubble constant H{t) = ^. The infinite universe'^ is 
treated by applying periodic boundary conditions to a 
finite cubic box of side L. 

We treat first the generic case that the N particles 
in the box are of equal mass and initially placed in any 
configuration which is a perturbed lattice, i.e., particles 
initially at the lattice sites of a perfect lattice subjected 
to some set of displacements. The right hand side of 
Eq. (1) may then be expanded perturbativcly in the rel- 
ative displacements of particles, about the perfect lat- 
tice configuration in which it is identically zero. Writing 
Xi(t) = R + u(R, i), where R is the lattice vector of 
the zth particle and u(R, t) its displacement from R, one 
obtains then, at linear order in this expansion, 

u{R,t) + 2Hu{R,t) = -^^I3(R-R')u(R',i). (2) 

" R' 

The matrix V is known in solid state physics, for any 
pair interaction, as the dynamical matrix (see e.g. [14]). 



1 The infinite sum for the force, as written in Eq. (1), is formally 
not well defined. It is implicit that it is regularized by the sub- 
traction of the effect of the mean density, the effect of which is 
taken into account in the expansion of the universe. In [13] we 
have studied the case of a static universe (i.e, the case where 
the scale factor a{t) = 1, in which the same regularization of 
this sum is applied) more extensively as it makes the comparison 
with the analogous condensed matter system much more trans- 
parent. The differences between the static and expanding case 
are, in any case, quite trivial for the PLT approximation. 



P,.(R^0) = Gm(|^-3:^) 



(3) 



where is the Kronecker delta. Note that a sum over 
the copies, associated with the periodic boundary condi- 
tions, is implicit in these expressions. 

The Bloch theorem for lattices tells us that V is di- 
agonalized by plane waves in reciprocal space. We can 
define the Fourier transform and its inverse by 



u(k,t) = ^e-*'"^u(R,t) 

R 

u(R,t) = l5]e^'^-^u(k,i) 



(5a) 

(5b) 



where the sum in Eq. (5b) is over the first Brillouin zone 
(FBZ) of the lattice, i.e., for a simple cubic lattice^ k = 
n(27r/L), where n is a vector of integers of which each 
component rii {i = 1, 2, 3) takes all integer values in the 
range < < N^'^ /2. 

Using these definitions in Eq. (2) we obtain 

u(k, t) + 2H{t)u{k, t) = - ^P(k)u(k, t) (6) 

where I?(k), the Fourier transform (FT) of 'D(R), is a 
symmetric 3x3 matrix for each k. 

Diagonalising 'D(k) one can determine, for each k, 

three orthonormal eigenvectors e„(k) and their eigenval- 
ues u;^(k) (n = 1, 2, 3). The latter obey a sum rule^ 



(7) 



where po is the mean mass density. 

Given the initial displacements and velocities at a time 
t = to, the dynamical evolution of the particle trajecto- 
ries is then given as 



u(R, t) = ^Yl l^^^^ + i)fi(k, to) 



„jkR 



(8) 



where the matrix elements of the "evolution operator" V 



^ The FBZ, for any cubic lattice, is defined as the set of N non- 
equivalent reciprocal lattice vectors closest to the origin k = 0. 

See [27] for the generalisation of the treatment described here to 
the case of a face centered cubic and body centered cubic lattice. 
^ This rule is known in the context of condensed matter physics 
[14] as the Kohn sum rule. 
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and Q are 



P^4k,i) =^[/„(k,f)(e„(k))^(e„(k)), (9a) 

n=l 

3 

Q^,(k,t) =^y„(k,t)(e„(k)V(e„(k)),. (9b) 



The functions [/„(k, t) and V^(k, i) are linearly indepen- 
dent solutions of the mode equations 



(10) 



chosen such that 



f/„(k,to) = 1, f/„(k,io) = 0, 
14(k,io) =0, •(4(k,to) = 1. 



(11) 



The determination of the evolution thus reduces to: 

• 1. diagonalization of the N 3x3 matrices ©(k), and 

• 2. solution of the 3N equations (10) and (11) for 

the mode functions. 

We note that 2?(k), and therefore the first step, depends 
only on the type of lattice, and not on the cosmology. The 
dependence on the latter [encoded in a{t)] comes only in 
the second step, in the mode functions Un{t) and Vn{t). 
Both steps are numerically straightforward (and require 
a number of operations proportional to A'^, quite feasible 
even with moderate computational power for as large 
as those used in the very largest current cosmological 
simulations). For certain cases, notably an EdS universe, 
the mode functions can be written analytically (see [13]). 
For completeness we give these expressions in App. A. 

Details of the diagonalization of 2?(k) are given in [13]. 
In Fig. 1 are shown results for a 16^ simple lattice. Specif- 
ically it shows the normalized eigenvalues 



(12) 



where po is the mean mass density (arbitrarily chosen at 
the time to). These are plotted as a function of the mod- 
ulus k of k, in units of the Nyquist frequency Un = j 
where t is the lattice spacing. The fact that the eigen- 
values at a given k do not have the same value is a direct 
result of the anisotropy of the lattice. Shown in the plot 
are also lines linking eigenvectors oriented in some chosen 
directions. This allows one to see the branch structure 
of the spectrum, which is familiar in the context of anal- 
ogous calculations in condensed matter physics (see e.g. 
[14]). The vectors k are those in the first Brillouin zone 
of the simple cubic lattice (as defined above). 

The number of particles in the simple cubic lattice con- 
sidered for the above calculation is small compared to 
that in current cosmological iV-body simulations, which 
now attain (e.g. [3]) values surpassing 1000^. While it 



11,0,0] 
[1,1,0] 

[1,1,1] 
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FIG. 1: Normalized eigenvalues (see text) for the eigenvectors 
of the displacement field on a simple cubic lattice with lattice 
spacing i. Each point corresponds to an eigenvector labelled 
by a vector k in the first Brillouin zone of the lattice, and fcjv 
is the Nyquist frequency. The different continuous lines link k 
vectors oriented in the indicated directions. Note that for the 
[1, O, O] and [1, 1, 1] directions the two tranverse eigenmodes 
are degenerate because of rotational symmetry. 



is quite feasible, as noted above, to do these calculations 
for such large particle numbers, it is of little interest to 
do so here: the effects we arc interested in depend essen- 
tially on the ratio of the wavenumber of a given mode to 
the Nyquist frequency. Increasing particle number only 
changes the density of reciprocal lattice vectors in these 
units (since 'k/kfq = rv/N^/^), just filling in more densely 
the plot of the eigenvalues in Fig. 1 but leaving its form 
essentially unchanged. 

In [13] we have investigated the domain of validity of 
PLT using numerical simulations, comparing the results 
of the evolution under full gravity with that obtained us- 
ing the formulae just given. We have considered "shuffled 
lattice" IC (in which a random uncorrclatcd displacement 
is applied to each particle of the lattice, corresponding to 
a power spectrum P{k) oc fc^ at small k), and cosmologi- 
cal type IC (a power spectrum cx k~'^). In both cases the 
full evolution is observed to be very well approximated 
until a time when the average relative displacements ap- 
proach the lattice spacing. 



B. Derivation of the fluid limit (Zeldovich 
Approximation) 

The continuum limit of this discrete analysis is easily 
obtained as follows. We consider the eigenmodes labelled 
by k. Sending the particle density to infinity, at fixed 
mean mass density po, we approach the asymptotic be- 
haviour seen in Fig. 1 as fc/fcjv 0. It is simple to show 
[13] that one has then one purely longitudinal mode (i.e. 
parallel to k) with e„(k) = 1, and two tranverse modes 
with e„(k) = 0. The evolution Eq. (8) may then be writ- 
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ten 



u(R,t) = u±(Ti,to)U^{t)+u^\(R,U,)U\\{t) (13) 



where 



u|, (R, to) = ^ 5](u(k, to) ■ k)k e'"'-^, (14a) 

k 

u±iK,to) = ^ ^[u(k,to) - (u(k,to) • k)k] e^"^-^, 

(14b) 



and analogously for the velocity fields, i.e., the longitudi- 
nal and transverse components of the initial displacement 
and velocity fields. The functions J7||(t) and V|| {t) are the 
linearly independent solutions of Eq. (10) for e„(k) = 1 
(i.e. w^(k) = -47rGpo) ), and Uj_{t) and V±{t) those 
for e„(k) = 1 (i.e. ti'^(k) = 0). We have observed in 
[12] that, for an EdS universe, this corresponds exactly 
to the solution found at linear order in a perturbative 
treatment of the equations for a self-gravitating fluid in 
the Lagrangian formalism in [15]. 

The solutions for the mode functions give, generically, 
a growing mode and a decaying mode. Since the for- 
mer always dominates after a sufficient time, there is an 
attractive asymptotic behaviour for the general solution 
Eq. (13) which may be written in the form 



u(R,i) = .9(t)q(R) 



(15) 



where g{t) is proportional to the purely growing mode, 
and q(R) is a longitudinal vector field. This is the Zel- 
dovich approximation (ZA) [28]. 



C. Evolution from ZA initial conditions 

In cosmological simulations the ZA is used to fix the 
IC, usually at a time at which the universe is well ap- 
proximated by an EdS cosmology. In this case we have 
(see App. A): 



Unit) 



3 /^y/' 

2 [to) ^ [to 



Vi\{t) = -to 



u±{t) 



2/3 



to 

V^{t) = 3to 



-1/3- 



The ZA may then be written 



3 / t 



u(R,i)= „ , , 

z \to 

v(R,i) =g(R,to)<, 



4/3 



g(R,io)io 



(16) 
(17) 
.(18) 



(19a) 

(19b) 



where g(R, t) and v(R, t) are the peculiar gravitational 
acceleration field and peculiar velocity field, respectively, 
defined by 



V = r — ax 
g = r — ax = r 

We therefore have 



U + 2-U 

a 



ux(R,io) = = vx(R,io) 
2 

V||(R,io) = — U||(R,io). 



(20) 
(21) 

(22a) 
(22b) 



Using these specific IC in the general expression for the 
PLT evolution given in Eq. (8), it is simple to show that 
this evolution may be written in the form 

= ^ E *)^^'"' ^o)^*'"" (23) 



where we have defined 



u(k, to) = k(k • u(k, to)) = k{t(k, to) , 



(24) 



using the fact that the initial displacement field is purely 
longitudinal. The set of vectors E(k, t) which encode the 
evolution in PLT of each mode of the displacement field 
are given by 

E(k,i) = V[C/„(k,t) + Ay„(k,t)](k • e„(k))e„(k) . 

(25) 

These can be calculated straightforwardly for the given 
lattice and cosmology. As emphasized above, the latter 
only enters in determining the mode functions, while the 
eigenvectors and eigenvalues are those of the simple cubic 
lattice widely used in cosmological simulations. 
The fiuid limit is expressed as 



E'*""'(k,t) = limE(k,t) 



mt) + —mt)]k (26) 



where, as discussed above, the limit is taken at fixed lat- 
tice spacing £, and C/||(i) and V||(<) are the solutions of 
Eqs. (10) and (11) for e„(k) = 1 (i.e. u)l{k) = -AwGpo). 
Thus we recover again the ZA, in which the initial dis- 
placements are simply amplified by the appropriate time 
dependent factor. 

From now on we will consider specifically the case only 
of the EdS cosmology. Since all the calculations we do 
with PLT are valid (and will be applied) only until the 
time of shell-crossing, this means we assume only that 
the cosmological model studied is well approximated by 
EdS from the starting red-shift until shell-crossing. This 
is almost always the case, notably in that of the currently 
favored ACDM model, as the starting red-shift is always 
well within the matter dominated era and the cosmolog- 
ical constant contributes to the expansion significantly 



6 



only at very low red-shift, well after shell crossing of the 
smallest included non-linear scales. For this case it is 
simple to verify, following the discussion in the previous 
subsection, that 



E""'^(k,t) 



a(t)k 



(27) 



where a = 1 at t = to- 

At sufficiently long times we note that the expression 
Eq. (25) will always be dominated by the most rapidly 
growing of the three modes at any given k, associated to 
the largest of three eigenvalues. Using the expressions in 
App. A for the EdS mode functions we can write 



E(k,i > to) 



1 2 + 3a+,,(k) 



3 (amax(k) + a,nax(k)) \ ^0 , 

X (k • e,„ax(k))e,„ax(k) (28) 



where the coefficients a^^xC^) are those given by the ex- 
pressions in Eqs. (A5) for the largest of the three eigen- 
values e„(k), with corresponding eigenvector e,„ax(k) ^. 



□ 

□ 

' o o 
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o 3" "<cos(e)<0.9 
□ 0.9<cos(e)<0.95 
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FIG. 2: Discreteness factor _Ddisp(k, t) quantifying the modi- 
fication with respect to the fluid limit (_Ddisp(k, i) = 1) of the 
power in the mode k of the evolved displacement field. The 
plot is given at at a = 5 (for a simulation starting at a = 1). 
See text for further details. 
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III. STATISTICAL MEASURES OF 
DISCRETENESS EFFECTS 

As we have explained, discreteness effects may be iden- 
tified, in the regime of validity of PLT, as the differences 
between the evolution described by PLT, and the fluid 
limit (FLT). For cosmological IC the effect of discrete- 
ness on any individual particle trajectory may thus be 
written as 



Audisc(R,i) - [E(k,t) 



E 



fluid/ 



ik R 



(k,t) {i(k,to)e"" 



(29) 

where the explicit expressions for the vectors E(k, t) and 
E''""^(k, i) have been given above. 

To quantify effects of discreteness in simulations we 
consider some statistical measures of the induced effects. 



A. The power spectrum of displacements 

The two point properties of the displacement field are 
specified by the matrix 



S^i,{'k,t) = w^(k, i)u*(k, t). 



(30) 



When the IC are set up, as just discussed, using the ZA, 
it follows from Eq. (23) that 

S^,{k,t) ^ E^{k,t)E,ik,t)\iiik,to)\'' (31) 



* We assume that k ■ emax(k) ^ 0, which is indeed the case for all 
k. 



5p^(k, to) = kf_,ku\u(k,to)\'' 



(32) 



where we used the fact that E(k, t) is real [and E(k, to) = 

!]■ ^ 

Let us consider the evolution of the trace of this matrix 
Poikjt) = TrS{k,t), for which we have that 



PoiKt) = |E(k,i)|2pc(k,io), 



(33) 



where, using the orthogonality of the eigenvectors, it fol- 
lows that 



3 r 



|E(k,t)p = ^ 



Un{k,t) + —Vn{k,t) 



(e„-k)2. (34) 



It is simple to verify [13] that the ensemble average 
of P£){k,t)/N over realizations of the IC is equal to 
the Fourier transform, as defined in Eq. (5), of the 
displacement-displacement correlation function (u(0, t) ■ 
u(R,t)) (where (..) denotes the ensemble average). 

To quantify the effects of discreteness, it is convenient 
to define the normalized quantity: 

i?.sp(k,t). 



Pl^"''^(fc,t) 

|E(k,t)|^ 

|Efluid(k,i)|2 

|E(k,t)|2 



a2(t) 



(35) 



where P|^""^(/c,t) is the trace of the power spectrum of 
displacements evolved in the fluid limit. 

In Fig. 2 is shown il'disp(k, t) for the value a = 5, i.e., 
after evolution from IC set at a = 1. Deviations from 
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unity arc a direct measure of the modification of the the- 
oretical (fluid) evolution introduced by the discreteness 
up to this time. Note that Z)disp(k, t) is plotted as a func- 
tion of k, each point corresponding to a different value 
of k. The fact that the evolution depends on the orien- 
tation of the vector k is a manifestation of the breaking 
of rotational invariance by the lattice discretisation. The 
three different symbols for the points correspond to three 
different intervals of the cosine of the minimum angle 9 
between the vector k and one of the axes of the lattice. 
We will return to this point in further detail below. 

We have chosen the value a = 5 because it is of the 
order of the typical one at which shell crossing occurs in 
A/^-body simulations, i.e., the initial amplitude in these 
simulations is typically such that this is the case^. As we 
have already noted, the discrete evolution and the fluid 
one arc described by diflferent exponents, so this figure of 
-Ddisp(k, a) will change as we change a. The resemblance 
between Fig. 2 and the optical branch of the dispersion 
relation in Fig. 1 allows us to infer simply this evolution: 
the reason for this resemblance is that already at the time 
chosen (a = 5), the expression (34) is well approximated 
by the regime t 3> io, given by Eq. (28), in which the 
most rapidly growing eigenmodes, which are those on 
the optical branches, dominate. Using Eq. (28) it is easy 
to verify that 



£)disp(k,t>to) ^ 
(k-emax(k))(2 



3amax(k)) 



3(Q!max(k) 



Clmax 



(k)) 



2[a-ax(k)-i] 



(36) 



For small k (compared to k^) we can simplify this ex- 
pression. In this case [13] the eigenvalues on the optical 
branch can be treated in a Taylor expansion about the 
fluid limit: 



e(k) « (1 - b(k)k 



(37) 



where 6(k) is a dimensionless coefficient of order unity 
which depends on the direction in reciprocal space. This 

expression is in fact [13], for the simple cubic lattice, a 



good approximation to the eigenvalues up to k 
Using it we have 



£'disp(k, a) 

and therefore 

£>disp(k,a > 1) ! 
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&(k)f2^2 

c(k))2a- 



X (k- 



(k. 



.(k))' 



|&(k)fc^ 



|b(k)fe^ 



k]M- 

(38) 
(39) 



^ This is true, for example, in the very widely used COSMICS 
package [29] for generating IC. The code in fact chooses the ini- 
tial red-shift given the number of particles, the comoving box-size 

and the cosmological model, using the criterion that the maximal 
mass density fluctuation of the realization on the grid is normal- 
ized to unity. This gives a variance at the lattice spa<;ing which 
is less than, but of the order of, unity. 




FIG. 3: Discreteness factor Z?disp(k, t) for a = 2, a = 5 and 
a = 10, averaged over bins of Ajkj = O.OSfcjv for a 64^ sim- 
ple cubic lattice. In each case is plotted also the equivalent 
average of the analytic approximation given in Eq. (38). 



In Fig. 3 is shown Ddisp(k, a), now averaged over all k in a 
narrow bin of fc, for different values of a. As anticipated 
we see that the discreteness effects, i.e., the deviation 
from unity of this quantity, grows as a function of time. 
Also shown are the corresponding averages of the approx- 
imate expression Eq. (38) ^\ We see that the agreement 
is very good, as expected, for suflaciently small k and a. 



B. The power spectrum of density fluctuations 

We have considered above the evolution of the two 
point properties of the displacement fields. Usually in 
the cosmological context one is more interested in a di- 
rect characterization of the evolution of the density fluc- 
tuations, notably through the power spectrum of density 
fluctuations or its Fourier transform, the reduced two- 
point (density-density) correlation function. In the cos- 
mological literature the relation between the two quan- 
tities is invariably derived in the continuous limit (see 
e.g. [30]). Instead we use here the more general result 
derived in [31], which gives the PS of density fluctuations 
of a discrete distribution of points generated by apply- 
ing a generic stochastic displacement field with two point 
correlation matrix 6'^i/(r) to a distribution with initial PS 
of density fluctuations Pj„(k): 

P(k) = ^ d^re 

- (2^)'^(k), 



(l + lm(r)) 
(40) 



^ We neglect here also, for simplicity, the time independent pre- 
factor. We will see below (cf. Eq. (55) and Fig. 7) that it is 



indeed very close to unity for k < k 
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where Cinij) is the Fourier transform of Pin(k), i.e., the 
reduced two point correlation function of the initial dis- 
tribution, and the integral is over all space (as the ex- 
pression applies in the infinite volume limit, i.e., to an 
infinite point distribution). 

Expanding the exponential factor e"'^'''^''!'^'*''^^)"'^'"'^'')! 
to linear order one obtains 



P(k) = Pi„(k)+fc^fc,5^,(k) (41) 
^3 I d3g5^,(q)[Pi„(k + q) - Pi„(k)] . 



Note that in reciprocal space we now have an integral, 
rather than a sum, as the density field is defined ev- 
erywhere in real space. For the case we are stiidying 
here, however, of a lattice with finite periodicity, the dis- 
placement field inherits the same periodicity and thus the 
function S^,y (k) is non-zero only on the reciprocal lattice. 
The integral therefore reduces again to a sum. 

In [5] we have studied in detail^ the domain of validity 
of this linearized approximation Eq. (41) to the full ex- 
pression Eq. (40). This depends in general on the shape 
of the input PS Pth(k), but for the PS in the range typical 
of cosmological models the criterion is just that the vari- 
ance of the relative displacement of particles be small 
compared to the distance separating them. This coin- 
cides precisely with the range of validity of PLT. Using 
then Eq. (41) with the PS of the displacements as given 
by Eq. (31), we obtain 



P(k, t) = Pi„(k) + (k • E(k, t)fPth{k) 



(42) 



+7^ / A(k-E(q,i))2:^[P,„(k + q)-Pi„(k)] 

where we have used the Poisson equation with Eqs. (19) 
to relate the input theoretical PS Pth{k) to the applied 
displacement fields through 



fc2|u(k,io)|2 = Pth(A;,io). 



(43) 



Eq. (42) is thus the expression for the evolved PS of 
the density fluctuations in and A^-body simulation (on a 
perturbed lattice) as given by PLT. The effects of dis- 
creteness described may be divided into two: 

• Effects encoded in the dependence of the function 
E(k, t) on the lattice spacing £, which give as we 
have discussed deviations from its fluid value ka(f). 
This is manifestly a dynamical effect of discretiza- 
tion on the lattice. 

• Effects encoded in the additional dependence on the 
lattice spacing £ of the first and third term on the 
right hand side of Eq. (42). This dependence comes 
through the initial PS Pj„(k). These terms describe 



the density fluctuations which are associated to the 
discrete sampling of the PS given by the second 
term on the right hand side of Eq. (42), on the 
lattice. They thus correspond to a static effect of 
the discretization. 

The first term on the right hand side of Eq. (42) de- 
pends only on the second effect — and is thus indepen- 
dent of time, describing simply the static contribution of 
the initial unperturbed lattice — while the second term 
depends only on the first effect. The third term, however, 
couples both effects. It describes, as we shall see, how 
purely discrete power, at large k is induced and evolves 
in time. 

Pi„(k), which is the PS of the lattice, is simply propor- 
tional to an infinite sum of delta functions at all points of 
the reciprocal lattice, i.e., at k = 2fcjvm where m is any 
vector of non-zero integers. This term thus vanishes in 
the first Brillouin zone [as defined after Eq. (12) above]. 
It is straightforward to show [5] that the same is true for 
the third (convolution) term in Eq. (42) , provided the in- 
put PS is appropriately set equal to zero outside the first 
Brillouin zone^. The second term on the right hand side 
of Eq. (42) is thus the only term which is non-zero inside 
the first Brillouin zone^. The non-trivial discreteness ef- 
fects described by the second and third terms on the right 
hand side of Eq. (42) thus contribute in non-overlapping 
regions of reciprocal space. 



1. PS inside first Brillouin zone 

As we have just seen the evolved PS inside the first 
Brillouin zone, to linear order in the input PS, is 

P(k,i) = (k-E(k,i))2Pth(fc). (44) 
To characterize the efi'ects of discreteness we thus define 



|k-E(k,i)|^ 



|k.Efluid(k,t)|2 

|k-E(k,t)|2 



(45) 



See, specifically, section IIIC and the appendices of [5]. 



* When applying the displacement field using Eq. (43) one can 
consistently include modes at larger wavenumbers. However this 
leads to aliasing efifects in the PS described precisely by this 
convolution term (i.e. to the appearance of unwanted power at 
long wavelengths). Sec [5] for further detail. 

^ As discussed in [5], this is not true for a glass prc-initial con- 
figuration. In this case the third term in Eq. (42) contributes 
at all k, and is proportional to at small k. See discussion in 
conclusions section. 
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FIG. 4: Discreteness factor Dsp{\i.,t) quantifying tlie modifi- 
cation with respect to the fluid hmit {Dgp(k,t) = f) of the 
power in the mode k of the evolved density fluctuation field, 
for k in the FBZ. The plot is given, as for Fig. 2, at a = 5 
(for a simulation starting at a = 1), for a 16^ simple cubic 
lattice. Also shown is a curve corresponding to the same quan- 
tity averaged over all k with k in one of bins of equal width 
Ak = O.OSfciv, for a 64"^ simple cubic lattice. (The points 
of this 64^ calculation, which are not shown, just trace more 
densely the behaviour of the points shown.) 



zone), and decreases as k does. At a = 5 the power in 
the largest mode is reduced to about one third of its fluid 
value, while around k = the fractional error varies 
from about -1-10% to more than —50%. At fc = fcAr/2 it 
varies from -1-5% to about —20%, while at A; = fcjv/4 the 
total spread is about 10%. 

In Fig. 4 is shown also an average of Dgp (k, a) (at 
a = 5) over narrow bins of equal width in k, for a larger 
64^ lattice. We see that this average is dominated, at 
this time, by the more numerous modes with growth co- 
efficients which are smaller than the fluid one. We see 
that this average describes at all k a net slowing down 
of the evolution of the power in the density fluctuations, 
ranging from slightly more than 30% at the Nyquist fre- 
quency, to 10% at half this frequency, and down to about 
3 — 4% at one quarter. It is straightforward to refine 
these estimates given the precise parameters of a simu- 
lation (i.e. the initial amplitudes to determine time of 
shell-crossing). In our conclusions we will discuss the 
importance of these effects, and how they might be cor- 
rected for in simulations. 



2. PS outside the first Brillouin zone 



where, using the definitions given above, we have 



|k.E(k,t)|2 



■ 3 

E 

.n=l 



3ii 



-K(k, i)](e„-k)2 



(46) 

instead of the expression in Eq. (34) for the analogous 
quantity for the PS of the displacement fields. 

In Fig. 4 is shown Dspl^,t), at a = 5 as for the anal- 
ogous Fig. 2 for the displacement fields in the previous 
subsection. The two figures are in fact very similar, par- 
ticularly at smaller k. The reason is the same one which 
explained the similarity between Fig. 2 and the optical 
branch of Fig. 1: the most rapidly growing mode on this 
branch already dominates at this time so that the differ- 
ence between the expression in Eqs. (34) and (46) reduces 
to a trivial time independent factor. Indeed we have now 



i:>5p(k,t > to) 



(2-f 3a 



+ 

max 



(k 
(k)) 



(k)) 

2 



(47) 



3(amax(k) 



(k)) 



which differs from Eqs. (36) only by the power of the 
product k-e„iax(k). We do not plot the analogous curves 
to those of Fig. 3 as the results look essentially the same. 

At a = 5 (i.e. at the time of shell crossing in a typi- 
cal cosmological A^-body simulation) our Fig. 4 is a plot 
of the fractional discrepancy between the theoretically 
evolved power (by FLT) and the power as evolved in the 
discretization of this system (by PLT). The fractional 
error introduced by the discretization is largest, unsur- 
prisingly, for the modes at the very largest wavenum- 
bers (fc = \/ikM, at the extremities of the first Brillouin 



For k outside the FBZ, we have in Eq. (42) only the 
non-trivial contribution from the third term. Taking an 
input theoretical PS of the form 



Ptn{k)^A¥'f{k/k,) 



(48) 



where f{k/kc) is a function which cuts off the PS at fc > 
kc, it is straightforward to show (see [5] for further detail) 
that this term can be written 



P!,'\k) = Ak'ky^Iik 



(49) 



where 



Ih-k| 



/(k) = ^[k.E(h-k,i)]' 

x/f^^T^ ) eFB2(h-k) (50) 



where h — 2fcArm, and the sum runs over all integer vec- 
tors m, and QpszO^—i^) is a three dimensional Heaviside 
function which is equal to unity inside the first Brillouin 
zone, and zero elsewhere. This cut-off is imposed, as dis- 
cussed above, in order to avoid aliasing effects. 

At a given k the sum /(k) can pick up contributions 
only from a single vector h, the one which lies inside the 
FBZ when translated by — k. It thus has the intrinsic 
periodicity of the lattice PS itself, i.e., /(k) ~ /(k + 
h), and it thus suffices to calculate it for vectors k such 
that kjsi < ki < 2k]\f. Just outside the FBZ it picks 
up contributions from h such that h — k lies just inside 
the FBZ, and thus of order /(fc^v/fcc), while at the larger 
values fc ~ 2fcjv it picks up contributions from h such that 



10 



0<k/k^<0.1 

0.1<k/k^<0.2 

0.2<k/kjj<0.3 



0.8 

cos(e) 



I 

0.9 



FIG. 5: Average value e of the normalized eigenvalues e„(k) = 
— ci;^(k)/47rGp(), averaged over the maximal eigenvalues at 
each k for k in the different ranges indicated, as a function of 
cos 6, the minimal angle between k and any axis of the lattice. 



h — k ^ 0, with an amplitude which depends strongly 
on the exponent n of the input power spectrum. For any 
n < 2, as is always the case in cosmological models, the 
amplitude diverges as one approaches k = 2fc„m (where 
Pin{k) is also divergent). 

The term in the PS given in Eq. (49) thus describes 
how the power at all wavenumbers in the input PS gives 
rise to power at small scales. Like the power Pi„(k) as- 
sociated with the unperturbed lattice, it is purely dis- 
crete. Differently from this latter term, which is time 
independent, it evolves in time as described by PLT. It 
thus describes evolving physical clustering in the simula- 
tions which is entirely an artifact of discreteness, rather 
than, as in the effects analyzed above, a modification due 
to discreteness of the non-trivial clustering of the fluid 
limit. Indeed it is interesting to note that, even if we 
take kc « fc^r, i.e., a cut-off in the input PS (as for 
example in typical in hot or warm dark matter simula- 
tions), this term is non-zero, and the growth can be well 
approximated by using the fluid limit for E(k, f) (since 
for the terms contributing in the sum |h — k|^ <C 1). We 
will discuss this term further in our conclusions below. 



C. Measures of anisotropy 

Let us consider in more detail the effects of anisotropy 
in the evolution induced by the lattice discretization. The 
different symbols indicated in Fig. 2 and Fig. 4 corre- 
spond, as has been noted, to the difFerent ranges of the 
cosine of the minimum angle 9 between the vector k and 
one of the lattice directions. This is, of course, just a 
direct result of the dependence of the eigenvalues on this 
quantity evident in Fig. 1, which is a manifestation of the 
anisotropy of the lattice. In Fig. 5 we plot the eigenvalue 
as a function of cos for the maximal eigenvalue (i.e. on 
the optical branch) averaged over different ranges of k. 



The degree of anisotropy over the range of orientations 
is roughly characterized by the slope of these curves (as 
zero slope corresponds to isotropy). We see clearly that 
the amplitude of the anisotropy thus decreases as the 
typical k in the chosen bin decreases, which is as ex- 
pected as we approach the fluid limit [for which s = 1]. 
We observe also a clear trend in the deviation from the 
fluid eigenvalue: this deviation typically increases as the 
orientation of the corresponding eigenvector goes further 
away from the axes of the lattice. This behaviour is, how- 
ever, not exactly followed, in particular around cos^ = 1 
. Indeed it is notable that the eigenmodes oriented par- 
allel, or very close to parallel, to the lattice planes have 
eigenvalues which arc larger than the fluid one. We note 
that the very largest eigenvalue in Fig. 1, associated with 
the largest values of Ddisp and Dsp in Figs. 2 and 4, cor- 
respond to an exactly longitudinal mode with k = kjsi 
and k parallel to the axes of the lattice. It thus de- 
scribes the motion of pairs of adjacent infinite planes to- 
wards one another The eigenvalue of this mode (sec 
Fig. 1) is e « 1.1, which gives [using Eqs. (48) and (A5)] 

There are clearly many different ways of quantifying 
these effects. One simple measure is the dispersion in 
the growth of the power in modes of the displacement or 
density fields. Considering the latter we define 



^Dsp{k,t) 



Dl,(k.l)-D.;{kA) 
D^\k,t) 



1/2 



|k-E(k,t)|2-|k.E(k,t)|2 
|k.E(k,t)|2 



(51) 

1/2 



where the second equality is valid in the FBZ. The aver- 
age, of any function X{\s.,t) on the reciprocal lattice, is 
defined as 



(52) 



k,|k|=fc 



Nk being the number of eigenmodes at a given k. 
ADsp{k,t) is evidently defined so that it is zero in the 
absence of anisotropy, and in particular in the fluid limit. 
Its value, averaged in narrow bins of k, is plotted in 
Fig. 6, for three different times corresponding to the 
scale-factors a indicated (where, as before, a = 1 corre- 
sponds to the beginning of the simulation). The results 
are very much as one would anticipate from Fig. 4: the 
anisotropy is greatest around k = fcjv where the spread 
in the growth factors is greatest, and the amplitude at 
any k grows as a function of a. Note that for a typical 
iV-body simulation, with a = 5 at shell crossing, the dis- 
persion in the growth factor at this time at A; = fcjv is of 
order 20%. 



Note that this is true for the case only of a lattice with N even. 
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a=2 

a=5 

a=10 



0.1 




FIG. 6: Discreteness factor ADsp{k,t) for a = 2, o = 5 and 
a — 10, averaged over bins of A|k| = 0. OS/cat for a, N = 64^ 
simple cubic lattice. 
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0.9 

0.8 I • • ' ' ' • • • 

0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 

FIG. 7: The asymptotic time independent anisotropy factor 

Ds,niso{k,t >> to) defined in the text averaged over bins of 
A|k| = O.OSfciv, for a 64'^ simple cubic lattice. 



One simple manifestation of the anisotropy introduced 
in the cvohition by the discretisation on a lattice was 
noted by Melott et al. in [19] : an initial perturbation 
described by a pure plane wave, which in the fluid limit 
should maintain its purely one dimensional character, can 
become three dimensional. The breaking of symmetry as- 
sociated is, up to numerical precision, entirely associated 
to the lattice. In [19] the authors studied this effect nu- 
merically, notably as function of the smoothing at small 
scales. One simple quantitative characterization of this 
effect is 



.^aniso(^5^) — 



|u(k,t) -k[k-u(k,t) 



(53) 



i.e., the fractional deviation of the displacement^^ direc- 
tion parallel to k. It follows from Eq. (23) that 



-^aniso(^; ^) — 



E^(k,t) 
(k-E(k,t))2 



(54) 



which is evidently unity in the fluid limit. For t ^ to it 
follows, from Eqs. (36) and (48), that 



-Daniso(k,t » to) 



1 



(k-e^ax(k))2 



(55) 



i.e., it is independent of time, depending only on the 

orientation of the eigenvector of the maximally growing 
mode with respect to k. This quantity, averaged over 
narrow bins of k, is shown in Fig. 7. We see the con- 
vergence to the fluid limit, in which emax(k) k, for 
fc/fcjv <C 1 and increasing deviations as k increases. 



In [19] the analogous quantity is considered, but for the velocity 
rather than the displacement. It is straightforward also to cal- 
culate this quantity; wc choose the displacement as it maJces the 
calculation a little simpler. 



IV. DISCRETENESS EFFECTS WITH A SMALL 
SCALE SMOOTHING 



In cosmological A^-body simulations a smoothing at 
small scales of the interparticle (Newtonian) potential is 
generally employed [1, 2]. "Small" in this context means 
that the characteristic scale for this modification of the 
potential is smaller or comparable to the average ini- 
tial inter-particle distance, i.e., approximately the lattice 
spacing £ in the case we are considering, of perturbations 
from a lattice. Almost always in fact this smoothing is 
taken to be considerably (i.e. one to two orders of magni- 
tude) smaller than £; an exception is in PM simulations 
in which the effective smoothing scale is normally equal 
to the scale £, as the potential is calculated on a mesh 
which typically is taken to coincide with the initial lattice 
configuration. There is in this case no correction to take 
into account forces from sub- mesh particles, as in P^M 
or tree codes (see e.g. [1, 2]). 

The motivation for the use of such a smoothing in this 
context is that it should reduce the effects of discreteness 
and make the simulation approximate better the theoret- 
ical (fluid-like) behaviour. The reason for this expecta- 
tion is that it reduces the impact of hard two-body col- 
lisions, which should be negligible in the fluid limit (in 
which the mean-fleld dominates). However, there is no 
demonstration (more than these approximate qualitative 
arguments) that such a smoothing has the effect of mak- 
ing a simulation a better approximation to the desired 
fluid limit. Indeed attempts in the literature to treat this 
question numerically (e.g. [16- 23]) lead to quite different 
conclusions. Notably the use of a smoothing smaller than 
£ has been contested by some authors [16-19], while such 
a practice is almost universal in large current A'^-body 
simulations. 

Let us consider what we can learn about this question 

using the methods developed in this paper. As was noted 
in Sect. II A, the perturbative analysis we have described 
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can be applied to any two-body interaction potential, and 
in particular to the case of a smoothed Newtonian poten- 
tial. We can simply repeat our calculation of the eigen- 
values and eigenmodes given the potential. The form of 
the smoothing employed varies and has evolved in time. 
Currently it is common to use a smoothing which leaves 
an exact 1/r potential for r > e, where e is the scale 
characteristic of the smoothing. If then e <C ^, as is typ- 
ically the case, this modification has no effect until two 
particles are closer than e. It therefore has no effect on 
our calculation (which breaks down when this condition 
is fulfilled anyway). Thus such a smoothing has no effect 
in moderating or otherwise the discreteness effects we 
have quantified. The fact that that this is the case illus- 
trates clearly an important point: the discreteness effects 
we have described are not two-body colUsional effects. 

In order for a smoothing of the potential to be of rel- 
evance to the effects we are describing it must modify 
the potential at the interparticle distance, as we do a 
perturbation expansion about the configuration with the 
particles at this distance. Let us consider therefore a sim- 
ple smoothing which does so, the so-called "Plummer" 
smoothing: 



(r) 



1 



(56) 



When e ~ £ it would be expected to, and indeed does, 
produce results similar to PM simulations (see e.g. [16, 
18]) 

For a given value of e we can calculate the spectrum 
of eigenvalues and eigenmodes. As one would anticipate, 
the modifications to them with respect to the case of 
pure gravity, are very small for e <C i!, and only becomes 
significant when we have e ~ f . In Fig. 8 we show the 
eigenvalues on the optical branch, for e = t/2. Also plot- 
ted, for comparison, is the same branch for pure gravity. 

The conclusions which follow from this figure are very 
clear. The effects of introducing a significant smoothing 
is to increase the average amplitude of discreteness ef- 
fects, while decreasing the effects of anisotropy. These! be- 
haviours are straightforward to understand qualitatively. 
The contribution to the total gravitational forces coming 
from nearby particles starts to be reduced significantly 
when the smoothing becomes of order the interparticle 
distance. Thus the growth factors at these scales de- 
crease relative to the case of full gravity, but they also 
show less anisotropy as the strong anisotropy in the dis- 
tribution of nearest neighbors is no longer felt as much. 

These effects can easily be further quantified by con- 
sidering the effects of varying the smoothing scale e on 
the quantities we defined above. In Fig. 9 we show the 
fluid-normalized discreteness factor .Ddisp [as defined in 
Eq. 45] both without any smoothing (as in Fig. 3) and 
for e = £/2. The result is self-evident given the results 
for the optical branch in Fig. 8 above: the smoothing 
increases the effects of discreteness on average. The ex- 
pected reduction in anisotropy as e increases is illustrated 
in Fig. 10, which shows the asymptotic value of the func- 




FIG. 8: Eigenvalues of the eigenmodes of the displacement 
fields for full gravity (top set), and with a smoothing at half 
the interparticle distance (see text), for a 64"* simple cubic 
lattice. For clarity only 0.5% of the eigenvalues are shown for 
fc/fejv > 0.2. 




FIG. 9: Discreteness factor Ddisp (k,t) for two different values 
of the smoothing parameter e (in units of I), for a = 2 and 
a = 10 for a 64'^ simple cubic lattice. An average over bins of 
width A|k| = O.OSfcAf has been performed. 



tion -Daniso [as defined in Eq. 53] for three values of e. As 
the latter parameter increases, we sec that on average the 
optical modes become more longitudinal which results in 
a reduced symmetry breaking effect in the propagation 
of plane waves. This is indeed the conclusion of the nu- 
merical study of this discreteness effect in [19] . We note, 
however, that the present analysis leads to a conclusion 
different to that of [19]: while the smoothing reduces this 
effect of discreteness, it actually increases the averaged 
effect of discreteness. By increasing the smoothing scale 
one obtains a better approximation to fluid behaviour 
only in the limited sense that one recovers better one 
qualitative feature of fluid behaviour, its isotropy. 
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FIG. 10: The asymptotic anisotropy factor £)aniso(fc,i 3> to) 
averaged over bins of A|k| = O.OSfejv, for three values of the 
smoothing parameter e, for a 64^ simple cubic lattice. 



V. PARAMETRIC AND LIMITING 
BEHAVIOURS 

We have given in the precedent sections various mea- 
sures of the modification of the fluid limit evohition en- 
gendered by discretization on a lattice. We have made 
use throughout of the result of [12, 13], which we have 
summarized in Sect. II B, that the fluid limit for the evo- 
lution of a given mode is recovered exactly when we send 
£ ^ while keeping k fixed. We have also seen that, 
keeping I fixed, the deviations of the evolution of this 
same mode from its fluid behaviour at shell crossing in- 
creases as the time of shell crossing docs. This implies 
that the evolution of cosmological simulations deviates 
completely from the desired theoretical evolution in the 
limit that the starting red-shift is increased keeping all 
other parameters fixed. Given these subtleties about how 
the continuous limit is recovered, it is instructive to con- 
sider these points in greater detail. 



A. The fluid limit 

Let us consider carefully in what limit the desired fiuid 
limit for the evolution of a cosmological iV-body simula- 
tion (in the class considered here) is recovered. 

Such a sinndation, periodic on a cube of side L, is a dis- 
cretisation of the continuum problem, which introduces, 
at least, the following parameters: 

• the lattice spacing £, given by i = L/N^''^, 

• the force smoothing scale e, and 

• the initial red-shift Zjnit- 

The interpretation of results of the simulation as physi- 
cal is justified if they approximate well those which would 

be obtained by simulations with suitable extrapolations 
of these parameters towards the theoretical continuum 



limit"'^^. The importance of understanding this limit pre- 
cisely is that it explains how this extrapolation of the 
parameters should be performed in practice. 

The input theoretical cosmological model, on the other 
hand, is characterized fully, for Gaussian initial fluctua- 
tions, by its PS. If the side L of the box is specified in 
physical units, the applicability of the ZA for setting up 
the IC then fixes, for given £, a minimal value of the ini- 
tial red-shift z-mit- This corresponds to an upper bound 
on the initial amplitude of the PS. The PS is, in gen- 
eral, a function which involves numerous physical scales 
characteristic of the cosmological model. It contains nec- 
essarily, however, at least one such scale: that charac- 
terising decay at large k of the PS. Indeed since the one 
point variance of the density fluctuations, which is pro- 
portional to the integral of the PS, is necessarily finite, 
one has that 



lim k^P{k) = . 



(57) 



There is thus always a physical scale given by a wavenum- 
ber kc, above which fluctuations of density decay faster 
than k^^. In simulations of the currently favored cold 
dark matter (CDM) models £ is always, because of nu- 
merical limitations, much larger than k~^ (i.e. kc is larger 
than the particle Nyquist frequency), and the displace- 
ment field generated is then cut-off in the first Brillouin 
zone rather than at kc (to avoid aliasing effects). In hot 
or warm dark matter (HDM or WDM) models, on the 
other hand, the (approximately exponential) cut-off is 
typically at a wavenumber considerably smaller than the 
Nyquist frequency. 

The exact fluid limit, in the regime of application of 
PLT, is that in which the corrections due to discreteness 
vanish, i.e., all the terms in Eq. (29) vanish exactly and 
the particles follow the exact FLT trajectories. It is clear 
that, strictly, this occurs only if we impose an abrupt 
cut-off kc in the input spectrum of modes, which remains 
fixed as the particle number increases, i.e., the exact fluid 
limit is kc£ ^ at fixed kc- We can, in practice, identify 
this kc with the (not necessarily abrupt) cut-off in the 
input PS: modes for k > kc do not contribute signiflcantly 
to the usually physical relevant quantities (such as mass 
variance) . 

To specify fully the limit we need to prescribe how the 
other parameters, ^;i„it and e, are treated. It is clear that 
it is sufficient to keep 2init fixed for what has been said 
above to hold. As for e, we have seen that it leads to 
a modification of the growth of modes with respect to 
their fluid evolution, so the latter will only be recovered 
if e ^ as £ — > 0, i.e., e — > in physical units. Given 



The results must of course also be, to a sufRcient approximation, 

independent of the box size L. Finite size effects of this type are 
not the subject of study here, as they are distinct from discrete- 
ness effects. Indeed such effects are equally present in the fluid 
limit (treated analogously with periodic boundary conditions). 
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that PLT is well defined for e = the order in which one 
takes the two limits is clearly not important. 

B. Long time/high initial red-shift limit 

We have seen that, because the modification of evo- 
lution associated to discreteness appears as a diff'erence 
in the exponents of growth of any given mode, the dis- 
crete system, with a given i, diverges from the fluid one 
when one considers the long time limit for any given k 
if I remains fixed. Indeed for k < fc^v we have scon, for 
example, that the statistical measures of discreteness ef- 
fects in the evolution of the PS of density fluctuations 
Dip is well approximated by Eq. (48). This expression 
diverges at large times arbitrarily far from the fluid limit, 
Dsp = 1. These PLT approximations remain good until 
close to shell crossing at the scale £ which occurs at a 
red-shift determined by the cosmological model (and in- 
dependent of the choice of the initial red-shift 2;init for the 
simulation). Thus, for a given physical scale, discreteness 
effects increase without limit as the starting red-shift of 
the simulation is increased, at fixed particle number. For 
asymptotically large values any physical quantity will be 
dominated by the modes with the largest exponents. In 
the case of the simple cubic lattice, as we have noted 
above, these modes correspond to the collapse of infinite 
planes along the axes of the lattice. 

It follows from this discussion that, in the fluid limit we 
deflned above, the order of the limits for N (or i) and for 
the duration t of the simulation (or 2;init)) is important, 
i.e., 

lim lim lim lim , (58) 

or, equivalently, 

lim lim lim lim . (59) 

t— »cxD Af— ►oo N—taot—>oc 

This result is in fact not surprising, at least for those fa- 
miliar with the use of the Vlasov limit to describe systems 
(e.g. plasmas) with long range interactions. It is in fact 
well known (e.g. [6, 7]) that there is a non-commutativity 
of the long time and large number limits in such systems. 
The Vlasov limit corresponds to taking the particle num- 
ber to infinity before the time. The long time limit of a 
system with fixed A'' is a different physical limit of the 
system's behaviour. Indeed a good example of this are 
precisely the two-body coUisional effects usually consid- 
ered as sources of discreteness error in A^-body simula- 
tion: these become relevant [32], in a uniform system 
with mean density po on a time-scale 1,2 ~ tdynN/ In N 
where tdyn = l/V^TrGpo- Thus clearly if one consid- 
ers the limit t — > oo at finite N one never recovers the 
Vlasov behaviour. Toy models in one dimension mani- 
festing this kind of behaviour have recently been studied 
in detail (see e.g. [33]). 

Practically this observation implies that at least one 
of the conditions for keeping discreteness effects under 



control in an A'^-body simulation is a constraint on the 
initial red-shift (i.e. given £, that zjnit be less than some 
value). We note that the z\-a\t is not a parameter which 
has been considered in discussions of discreteness effects 
in N-hody simulations in the literature (e.g. [16-19, 23]). 

C. N dependence of effects 

It is instructive also to consider more explicitly the N 
dependence of the effects we have described. A charac- 
teristic time scale for these effects, somewhat analogous 
to the one just refered to for two-body relaxation, can be 
derived by considering the evolution we have described 
here in static space (i.e. with a(t) = 1). The only modifi- 
cation with respect to the EdS case we have treated here 
is in the mode functions, which become simply growing 
and decaying exponentials. Fluctuations on all but the 
shortest scales have eigenvalues which may be well ap- 
proximated by Eq. (37), or 

.;(k) [1 - ^fc^£^)]C (60) 

where a'(k) = f^y^ is the fluid behaviour [and 6(k) is a 

constant of order unity] . For any quantity dominated by 
modes around k the deviation from fluid behaviour will 
become significant when 

The characteristic time scale for the discreteness effects 
to dominate may thus be estimated as 

UUk) ~ ^ oc (62) 

where the last proportionality is inferred for a given k 
(fixed in units of the box size as the particle number N 
increases). Note that here N is the number of particles 
sampling the scale ^ rather than the total number 
of particles in the system. Thus the time scale for the 
discreteness effects to dominate is a function of the scale 
considered rather than a global timescale for the system. 

VI. DISCUSSION AND CONCLUSIONS 

Wc have described in this paper the application of 
a formalism developed in [12, 13] to determine the er- 
rors in Af-body simulations at early times, i.e., up to 
shell crossing when the formalism breaks down. We have 
given an explicit expression, Eq. (29), for the error in 
the trajectory of any particle, assuming only that IC 
arc set up in the standard way using the ZA. We have 
then defined, evaluated numerically for a typical red-shift 
of shell- crossing, and given approximate analytic expres- 
sions for, fimctions which characterize these discreteness 
errors for any given mode of the displacement or density 
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field. Wc have also discussed the cfi'ects of a smooth- 
ing of the force at small scale on these effects, as well as 
the precise limits in which the fluid limit of the system's 
evolution is recovered. 

We now discuss both the immediate practical relevance 
of these results — their application in quantitifying a 
source of systematic error in simulations — as well as 
their relevance to the broader question of discreteness 
error in iV-body simulation, for other kinds of initial con- 
figurations and beyond shell crossing. 



A. Application of results 

We have seen that, at a red-shift of shell crossing typi- 
cal for a cosmological simulation (with IC using e.g. the 
code of [29]), errors larger than a few percent have accu- 
mulated for all wavenumbers larger than about a quar- 
ter of the Nyquist frequency. Such modes are usually 
included in the IC in simulations of CDM type mod- 
els, which typically cut the input PS at the edge of the 
first Brillouin zone. Indeed this means that for most of 
the modes included errors are larger : the modes with 
k > k^/A represent more than 99% of the sampled modes 
while those with k > k^ f'^, which have errors exceed- 
ing 25%, about 93%. Any physical quantities at the end 
of the simulation which depends on the power in these 
modes must inherit at least this error. Such power is in- 
cluded usually as one is indeed interested in such quan- 
tities. 

What precisely the induced error is will depend of 
course on the quantity considered, and we do not attempt 
here to describe such cfi'ects exhaustively. The results 
given in this paper allow one to calculate them, and thus 
in principle correct for them, in a very straightforward 
manner: it suffices to calculate the function £)(iisp(k, apu) 
as described in Sect. Ill A. We have introduced here the 
scale factor labelled Opit to indicate that up to which 
PLT is a good approximation. As we have discussed, 
and shown in detail by numerical study in [13] , this cor- 
responds approximately to shell crossing. Its precise def- 
inition is of course somewhat subjective as it involves 
specifying how strong a deviation defines it (and in what 
quantity the deviation is measured). Here it can be taken 
either as a theoretically calculated shell crossing at the 
scale £ or, better, as determined a posteriori from the sim- 
ulation (e.g. by matching the evolution of the variance 
of displacements to that described by PLT, as reported 
in [13]). 

The effects of discreteness up to shell crossing may thus 
be considered as equivalent to a k-dependent misrepre- 
sentation of the input displacement field (at a = 1), by 
a factor in amplitude of 1/ \/ £'ciisp(k, flpit). The effective 



initial amplitude is, for almost all modes, somewhat less 
than the theoretical amplitude. This could, in principle, 
be "corrected" for by multiplying each mode k of the in- 
put displacement field, at the chosen Zjnit (corresponding 
to a = 1), by the factor l/y^ £'disp(k, Opit). Such a correc- 
tion would ensure that the full linearized evolution PLT 
coincides, at a = Opit, with the analogous (linearized) 
fluid result. 

It is not clear, however, that such a correction of the 
initial amplitude will in fact produce a more accurate re- 
sult than simply starting the simulation at a = flpit, using 
the ZA to set up the IC at this time (and possibly im- 
proved taking higher order corrections in Lagrangian per- 
turbation theory into account). Applying the discrete- 
ness correction to the IC at a = 1 (< apu) means that 
the linear evolution is offset from its theoretical evolution 
for a < ttpit- It is possible that the integrated non-linear 
effect of this offset will be greater than the error involved 
in using the ZA directly at a = Opit when the linear ap- 
proximation it is based on is breaking down. In principle 
there exists a starting red-shift which is optimal in the 
sense that it minimizes the combined errors due to non- 
linearity (which decreases as Zinit increases^^) and errors 
due to discreteness (which increases, as we have seen, as 
-Zinit does). We note in this respect that the perturbative 
formalism used here, which described the fully discrete 
problem, can be extended beyond the linear order we 
have used. Such a treatment should allow a determina- 
tion of this optimal red-shift. 

These comments apply to IC as they arc applied in 
CDM type simulations. In HDM or WDM simulations 
the physical cut-off in the initial PS is, as has been noted 
above, typically such that kc « Aijv- Thus the modes 
of the displacement field which have the greatest mod- 
ification of their evolution due to discreteness arc not 
present. For modest typical values of Zinit the effects 
just discussed are therefore negligible. Above, however, 
we discussed only the effects of the modification of the 
dynamical fluid evolution, and not the additional effects 
which have been described in Sect. IIIB, in the PS of 
density fluctuations outside the FBZ. We have shown in 
our discussion of these terms that this additional, purely 
discrete, contribution in the PS is present even if the in- 
put PS is limited to a region well inside the FBZ, as in 
HDM or WDM simulations. In the approximation that 
the small k modes grow according to fluid theory (i.e. are 
linearly amplified in proportion to the scale factor a in 
the EdS model) it is simply an extra contribution in the 
initial PS which is also linearly amplified. It describes 
the growth of power at small scales due to the coupling 
between long wavelength power and that intrinsic in the 
initial lattice. Wc note that this physical mechanism gen- 
erating "artificial" discrete power at scales of order the 



The fraction of modes with k < fco is, approximately, the volume For evaluation of the errors in the fluid limit due to non-lineajrity, 

of a sphere of radius feo divided by the volume of the FBZ. see [8, 9]. 
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initial intcrparticlc distance has been observed and stud- 
ied in numerical simulations of HDM/ WDM [25, 26]. Our 
treatment gives an analytic description of this process (up 
to shell crossing), for the case of an initial simple cubic 
lattice. 

In this latter context comparative studies have been 
made of this initial configuration with ones generated 
starting from the "glass" configurations sometimes used 
as an alternative [34]. The formalism we have used in 
this paper can in principle be generalized to treat such a 
starting configuration, and thus to compare analytically 
the discreteness effects up to shell crossing in the different 
cases. There is, however, a considerable technical compli- 
cation: the Bloch theorem does not apply and the eigen- 
vectors of the displacement field about the (now approxi- 
mate) zero force configuration arc no longer plane waves. 
Thus the diagonalization problem cannot be simplified 
in the same way, and a fully numerical solution would 
be required to determine the 3N eigenmodes. We expect 
to find only small quantitative differences: the origin of 
the discreteness effects we have quantified on a lattice 
arise essentially, as we have discussed just above, from 
the sampling of the continuous density field by particles. 
The magnitude of the effects will depend essentially on 
the average particle density, and not on the details of the 
arrangement of the points in the distribution. 

We note, on the other hand, that the treatment of the 
problem on other cubic lattices is straightforward. In 
[27] the cases of a body-centered cubic and face-centered 
cubic lattice are solved. We will exploit this treatment 
of discreteness effects on a range of different "pre-initial" 
lattices in a forthcoming numerical study of discreteness 
effects in the non- linear regime [35]. 



B. Discreteness in the non-linear regime 

The perturbative treatment presented here can, as has 
been noted, be extended to beyond linear order. While 
such an approach may lead to useful insights into the 
interplay of discreteness effects and non-linearity, it can 
extend at most as far as shell crossing. The current for- 
malism therefore allows us to say nothing quantitatively 
about discreteness effects in the fully non-linear regime. 
It can however, as we now discuss, give us some qualita- 
tive insight which may be useful in attempts to under- 
stand and quantify such effects. 

Discreteness effects, we recall, are all those effects 
which lead to differences between the evolution of sta- 
tistical quantities in A''-body simulations and that in the 
theoretical Vlasov-Poisson limit. We have underlined 
that the effects studied here, which are manifestly of 
this type, are different in nature to those usually empha- 
sized in this context: (i) two-body coUisional effects (sec 
e.g [21, 22, 24]), and (ii) the limitations imposed by the 
particle discretisation on the representation of the initial 
theoretical PS (described analytically in [5]). The mod- 
ifications of the continuum evolution we have described 



here depend explicitly on the dimensionless ratio k£, and 
decrease in amplitude as ki decreases. They can thus 
be well described, and distinguished from these other ef- 
fects, as dynamical sparse sampling effects, i.e., they are 
modifications of the evolution of fluctuations which arise 
from the fact that a given scale is sampled more sparsely 
by the "macro-particles" of an iV-body simulation than 
by the microscopic particles of a fluid. 

Clearly one would expect that such a physical ef- 
fect should be present also in the evolution after shell- 
crossing, in addition to two-body collisional effects and 
any effects inherited from the IC. The difficulty in as- 
sessing its importance, at a given scale and time, is that 
it is not clear what parameter should play the role of i 
after shell crossing. The results of simulations are usu- 
ally interpreted as physical down to scales ~ e, where 
e <C ^. This is typically of the order of the interparticle 
distance in the densest regions of the final configurations 
(at the centers of halos). This extrapolation would ap- 
pear then to assume, roughly, that the role of d in the 
non-linear regime is played by the minimal interparticle 
distance in the evolved simulation. This seems, given the 
physical nature of the effect, an extremely optimistic as- 
sumption. One would expect fluctuations at signiflcantly 
larger scales in any less dense region (i.e. most of the 
volume) to be subject to much greater discreteness error 
of this type. In very low density regions, notably, it is 
clear that all but the very longest wavelengths (compared 
to the locally very large interparticle separation) will be 
very poorly approximated. 

Another related point which emerges from our analysis 
of the evolution up to shell crossing is the subtlety in- 
volved in deflning the correct continuum limit. We saw, 
in particular, that an extrapolation to high initial red- 
shift at constant particle number (and force smoothing) 
leads to a divergence from the desired limit. We note 
that this is an unexpected behaviour: variation of this 
parameter has been considered as relevant only to the 
accuracy with which the IC are represented [8, 9], as 
the Zcldovich approximation becomes exact only in the 
limit that 2;init — > oo. Our results indicate therefore that 
numerical tests which vary A'' to test for discreteness ef- 
fects should best be done at fixed zjnii. Further we have 
seen that the continuum limit is attained strictly only 
in the limit that fcc^ — > where kc is the cut-off in the 
input PS. While we have been able to calculate correc- 
tions outside this regime (i.e., for CDM type simulations 
which, as noted above, are performed with kj. ^ 1) it 
is not evident that it will be possible to do so in the 
fully non- linear regime. This suggests that, to calculate 
robust discreteness errors in the latter regime, it may in- 
deed be necessary to extrapolate numerically to a regime 
in which £ (as argued notably in [16-19]) — or what- 
ever length scale plays its role in the non-linear regime 
— is much less than all physical scales at which cluster- 
ing is to be determined. A related discussion of some 
of these issues, based on numerical study of simplified 
"cosmological-like" simulations, can be found in [36, 37]. 
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Wc will explore them further for the specific question 
of quantification of errors in this regime in forthcoming 
publications. 
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APPENDIX A: EIGENMODES IN AN EDS 
UNIVERSE 

The evolution of the displacement field depends, as 
we have emphasized, on the cosmological model only 
through the mode functions [/'„(k, and V^(k, t), which 
are solutions of Eq. (10) satisfying the conditions in 
Eq .(11). In the the case of an EdS universe the former 
becomes 

/(k, t) + |/(k, t) = Ae4k)/„(k, i), (Al) 

where, as above, e„(k) = — a;^(k)/47rGpo, and we have 
used 
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It follows that 
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an(k)= g [Vl + 24e„(k)-l 



a+(k) = - [Vl + 24e„(k) + 1 



(A5a) 

(A5b) 



Note that in the fluid limit we obtain a_ = 2/3 and 
q:+ = 1 for the longitudinal mode (e„ = 1) and a_ = 
and a+ = 1/3 for the tranverse modes (e„ = 0). If 
e„(k) > the solution arc a power- law growing mode 
and a power-law decaying mode. If > e„(k) > —1/24, 
there are two decaying modes. Finally, if e„(k) < —1/24, 
the solution is oscillatory and can be written as 
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